home *** CD-ROM | disk | FTP | other *** search
/ BBS in a Box 3 / BBS in a box - Trilogy III.iso / Files / Prog / U-Z / VideoToolBox Folder / Utilities / Quick3 / Weibull.c < prev   
Encoding:
C/C++ Source or Header  |  1993-02-23  |  1.5 KB  |  45 lines  |  [TEXT/KAHL]
  1. /*
  2. Weibull.c
  3. Copyright (c) 1990, 1991, 1992 Denis G. Pelli 
  4. A psychometric function for use by PsychometricFit() and Quick3.
  5. Computes a Weibull function. E.g.
  6. alpha=0.01; threshold
  7. beta=3.5;    steepness
  8. gamma=0.5;    guessing rate
  9. delta=0.0;    rate of random guesses at high, suprathreshold contrasts
  10. P=Weibull(contrast,&myParamRecord);
  11. Note: Weibull's prototype must be consistent with PsychometricPtr in Quick3.h
  12.  
  13. HISTORY:
  14. Fall'89        dgp    wrote it, but didn't use it.
  15. 4/8/90        dgp    polished it. Added ILLEGAL_PARAMETERS return value to deal gracefully
  16.                 with the fact that general purpose minimization routines will quite 
  17.                 naturally try illegal values of the parameters and that the appropriate 
  18.                 response is that such parameters produce a very bad fit, inducing the 
  19.                 minimization program to stay within legal values of the parameters.
  20. 10/29/90    dgp    Tidied up comments.
  21. 3/11/92        dgp For speed I changed pow(10,...) to exp10(...), which is defined in
  22.                 Quick3.h.
  23. */
  24. #include "Quick3.h"
  25.  
  26. double Weibull(double contrast,paramRecord *paramPtr)
  27. {
  28.     double p,exponent,alpha,beta,gamma,delta;
  29.  
  30.     alpha = exp10(paramPtr->logAlpha);
  31.     beta=paramPtr->beta;
  32.     gamma=paramPtr->gamma;
  33.     delta=paramPtr->delta;
  34.     if(beta<0.0 || gamma<0.0 || gamma>1.0 ||
  35.         delta<0.0 || delta>1.0) return ILLEGAL_PARAMETERS;
  36.     contrast /= alpha;
  37.     exponent = pow(contrast,beta);
  38.     if ( exponent > 30.0 )
  39.         p = 1.0;    /* This avoids floating underflow */
  40.     else
  41.         p = 1.0 - (1.0 - gamma) * exp(-exponent);
  42.     p = (1.0-delta)*p + delta*gamma; 
  43.     return p;
  44. }
  45.